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EXPLICIT TIME-STEPPING FOR STIFF ODES* 

KENNETH ERIKSSON 1 " , CLAES JOHNSON"!", AND ANDERS LOGG* 

Abstract. We present a new strategy for solving stiff ODEs with explicit methods. By adap- 
tively taking a small number of stabilizing small explicit time steps when necessary, a stiff ODE 
system can be stabilized enough to allow for time steps much larger than what is indicated by classi- 
cal stability analysis. For many stiff problems the cost of the stabilizing small time steps is small, so 
the improvement is large. We illustrate the technique on a number of well-known stiff test problems. 
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1. Introduction. The classical wisdom developed in the 1950s regarding stiff 
ODEs is that efficient integration requires implicit (A-stable) methods, at least outside 
of transients, where the time steps may be chosen large from an accuracy point of 
view. Using an explicit method (with a bounded stability region) the time steps have 
to be small at all times for stability reasons, and thus the advantage of a low cost per 
time step may be counterbalanced by the necessity of taking a large number of steps. 
As a result, the overall efficiency of an explicit method for a stiff ODE may be small. 

The question now is if it is possible to combine the low cost per time step of an 
explicit method with the possibility of choosing large time steps outside of transients. 
To reach this favorable combination, some kind of stabilization of the explicit method 
seems to be needed, and the basic question is then if the stabilization can be realized 
at a low cost. 

The stabilization technique proposed in this note relies on the inherent property of 
the stiff problem itself, which is the rapid damping of high frequencies. This is achieved 
by stabilizing the system with a couple of small stabilizing (explicit Euler) steps. We 
test this idea in adaptive form, where the size and number of the small time steps 
are adaptively chosen according to the size of different residuals. In particular, we 
compute rates of divergence to determine the current mode A of largest amplification 
and to determine a corresponding small time step k « 4 with high damping. We 
test the corresponding adaptive method in the setting of the cG(l) method with 
fixed point iteration, effectively corresponding to an explicit method if the number of 
iterations is kept small. We show in a sequence of test problems that the proposed 
adaptive method yields a significant reduction in work in comparison to a standard 
implementation of an explicit method, with a typical speedup factor of ~ 100. We 
conclude that, for many stiff problems, we may efficiently use an explicit method, if 
only the explicit method is adaptively stabilized with a relatively small number of 
small time steps, and so we reach the desired combination of a low cost per time step 
and the possibility of taking large time steps outside transients. It remains to be seen 
if the proposed method can also be efficient in comparison to implicit methods. 

We have been led to this question in our work on multi-adaptive cG(q) or dG(g) 
ODE-solvers based on Galerkin's method with continuous or discontinuous polyno- 
mials of degree q, where individual time steps are used for different components (see 
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[SJIH])- These methods are implicit, and, to realize efficient implementations, we need 
to use fixed point iteration with simple preconditioners. With a limited (small) num- 
ber of iterations, these iterative solvers correspond to explicit time-stepping, and the 
same question of the cost of stabilization arises. 

The motivation for this work comes also from situations where for some reason 
we are forced to using an explicit method, such as in simulations of very large systems 
of chemical reactions or molecular dynamics, where the formation of Jacobians and 
the solution of the linear system become very expensive. 

Possibly, similar techniques may be used also to stabilize multigrid smoothers. 

2. Basic strategy. We consider first the test equation: Find u : [0,oo) — > M. 
such that 

u(t) + Xu(t) =0 for t > 0, 

(2.1) W W 

V ' u(0) = u°, 

where A > and u° is a given initial condition. The solution is given by u(t) — 
exp(— Xt)u°. We define the transient as {t > : \t < C} with C a moderate constant. 
Outside the transient, u(t) = exp(—Xt)u° will be small. 
The explicit Euler method for the test equation reads 

U n = -knXU" 1 - 1 + U 71 - 1 = (1 - k n \)U n -\ 

This method is conditionally stable and requires that k n X < 2, which, outside of 
transients, is too restrictive for large A. 

Now let K be a large time step satisfying KX > 2 and let A; be a small time step 
chosen so that kX < 2. Consider the method 

(2.2) U n = (l-kX) m (l-KX)U n -\ 

corresponding to one explicit Euler step with large time step K and m explicit Euler 
steps with small time steps fc, where m is a positive integer to be determined. Alto- 
gether this corresponds to a time step of size k n = K + mk. Defining the polynomial 
function P(x) = (1 — 9x) m (l — x), where 9 = j^, we can write method (|2.2[) in the 
form 

U n = PiKX)^- 1 . 

We now seek to choose m so that |P(ifA)| < 1, which is needed for stability. We thus 
need to satisfy 

|l-fcA| m (^A-l)<l, 

that is, 

log(ia-l) \og(KX) 



(2.3) m > 



log|l-fcA| 



with c = fcAa moderate constant; for dcfinitcncss, let c — 1/2. 

We conclude that m will be of moderate size, and consequently only a small 
fraction of the total time interval will be spent on time-stepping with the small time 
step k. To see this, define the cost as 
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that is, the number of time steps per unit time interval. Classical stability analysis 
with m = gives 

(2.4) a = l/k n = A/2, 

with a maximum time step k n = K = 2/ A. Using (|2.3[) we instead find 

for if A ^> 1. The cost is thus decreased by the cost reduction factor 

21og( J ftTA) log(iTA) 
~cK\ KX ' 

which can be quite significant for large values of KX. 

A similar analysis applies to the system it + Au — if the eigenvalues {Xi}fLi °f 
A are separated with < Ai < • • ■ < Aj-i < 2/K and 2/K <C Aj < • • • < A at. In this 
case, the cost is decreased by a factor of 

(26) 21og(^A 4 ) log(^A 4 ) 



ci^A, XA, 

In recent independent work by Gear and Kevrekidis [2] , a similar idea of combining 
small and large time steps for a class of stiff problems with a clear separation of slow 
and fast time scales is developed. That work, however, is not focused on adaptivity 
to the same extent as ours, which does not require any a priori information about the 
nature of the stiffness (for example, the distribution of eigenvalues). 

3. Parabolic problems. For a parabolic system, 

u(t) + Au(t) = for t > 0, 
(3.1) W W . 

u(0)=u°, 

with A a symmetric positive semidefinite N x N matrix and u° € a given initial 
condition, the basic damping strategy outlined in the previous section may fail to 
be efficient. This can be seen directly from (|2.6j) : for efficiency we need KXi ^> 2, 
but with the eigenvalues of A distributed over the interval [0, Aat], for example with 
Ai <~ i 2 as for a typical (one-dimensional) parabolic problem, one cannot have both 
Aj_i < 2/K and KX, > 2! 

We conclude that, for a parabolic problem, the sequence of time steps, or equiva- 
lently the polynomial P(x), has to be chosen differently. We thus seek a more general 
sequence k\, ■ ■ ■ , k m of time steps such that |P(x)| < 1 for x € [0, ifAjy], with 

P ( X ) = ( 1 -t^.)...(i kmX 



K \ K 



and K a given maximum step size. 



3.1. Chebyshev damping. A good candidate for P is given by the shifted 
Chebyshev polynomial of degree m, 
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This gives P c (0) = 1 and |P c (x)| < 1 on [0,KXn] (see Figured]). A similar approach 
is taken in [llj . 

Analyzing the cost as before, we have a = m/ (ki + • • ■ + k m ), with 

h = 2/(Ajv(1 - s m+ i_j)), 

and Si the ith zero of the Chebyshev polynomial T m = T m (s), given by Sj = cos((2« — 
l)7r/(2m)), i = l, . . . ,m. It follows that 

mX N /2 m\ N /2 

a = — — = - — = An Am. 

1/(1- «!) + ■■• + 1/(1 -s m ) m 2 N/ 

The reduction in cost compared to Xn /2 is thus a factor l/m. A restriction on the 
maximum value of m is given by k m = 2/(Aat(1 — si)) < K\ that is, 

2 2 

K > ; ; rr-T ~ TTV, ?TT = 16m, 2 /( A MTT 2 ) 

- Ajv(l - cos(7r/(2m))) AAr7r 2 /(8m 2 ) /yN ' 



for m ^ 1. With m = {-k/A)^/K\n , the cost reduction factor for Chebyshev damping 
is thus given by 



l/m = — — (K\ 



N 



-1/2 



3.2. Dyadic damping. Another approach is a modified version of the simple 
approach of small and large time steps discussed in the previous section. As discussed 
above, the problem with this basic method is that the damping is inefficient for 
intermediate eigenmodes. To solve this problem, we modify the simple method of 
large and small time steps to include also time steps of intermediate size. A related 
approach is taken in [3]. Starting with a couple of small time steps of size k, we thus 
increase the time step gradually by (say) a factor of 2, until we reach the largest step 
size K. Correspondingly, we decrease the number of time steps at each level by a 
factor of 2. The polynomial P is now given by 

*«-.n(i-&)n(i-£' J 

J=q+1 i=Q 

where p and q < p are two integers to be determined. We thus obtain a sequence of 
increasing time steps, starting with 2 q time steps of size k = 1/Ajv, 2 9 ~ 1 time steps of 
size 2k, and so on, until we reach level q where we take one time step of size 2 q k. We 
then continue to increase the size of the time step, with only one time step at each 
level, until we reach a time step of size 2 p k — K . 

With p given, we determine the minimal value of q for |Pd(a;)| < 1 on [0, K Aat] 
for p = 0, 1, . . . , and find q{p) ~ |(p — 1); see Table [T] The cost is now given by 

(2' + — + l) + (p-«/) Xn^ 1 -l) + (p-q)) 



j^(2i(q + 1) + (29+ 1 + • • • + 2Pj) 2i(q + 1) + (2P+ 1 - 21+ 1 ) 

Ajv/2 ^V -p+2(p-l)/3+l _ A 



AN_ 2 -p+2(p-l)/3+l _ A N 2 -(p-i)/ 3 



2P-9-1 2 2 
Since p = log{K Ajv)/log(2), the reduction in cost for dyadic damping is then a factor 

1 



(A-Aat/2)V3 

which is competitive with the optimal cost reduction of Chebyshev damping. 
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Table 1 

Number of required levels with multiple zeros, q, as a function of the total number of levels, p. 



p 


1 


2 


3 4 


5 


6 7 8 


9 


10 


11 12 


13 14 


15 


16 


1 








1 2 


3 


3 4 4 


5 


6 


7 8 


8 9 


10 


10 



4. Comparison of methods. We summarize our findings so far as follows. The 
reduction in cost for the simple approach of large unstable time steps followed by small 
stabilizing time steps (outlined in section [2]) is a factor log(ATAjv) / '(KXn), and thus 
this simple approach can be much more efficient than both of the two approaches, 
Chebyshev damping and dyadic damping, discussed in the previous section. This 
however requires that the problem in question has a gap in its eigenvalue spectrum; 
that is, we must have a clear separation into small (stable) eigenvalues and large 
(unstable) eigenvalues. 

In the case of a parabolic problem, without a gap in the eigenvalue spectrum, 
the simple approach of section [5] will not work. In this case, Chebyshev damping and 
dyadic damping give a reduction in cost which is a factor (KXn)- 1 / 2 or (KXn)- 1 / 3 , 
respectively. The efficiency for a parabolic problem is thus comparable for the two 
methods. Since the dyadic damping method is a slightly modified version of the simple 
approach of section consisting of gradually and carefully increasing the time step 
size after stabilization with the smallest time step, thus being simple and flexible, we 
believe this method to be advantageous over the Chebyshev method. 

As a comparison, we plot in Figure [TJ the shifted Chebyshev polynomial P c (x) 
and the polynomial Pd(x) of dyadic damping for KX N = 64. With KXn = 64, we 
obtain m = (ir/ 4)\/KXn ~ 6 for the Chebyshev polynomial. For the dyadic damping 
polynomial, we have p — log(-ftfAAr)/log(2) = 6, and from Table [TJ we obtain q = 3. 

As another comparison, we plot in Figure [5] the stability regions for the two 
polynomials P c (z) and Pd(z) for z G C. In this context, the two polynomials are 
given by 

(4.i, *w-n(i + i^) 

for Si = cos((2« — l)7r/(2m)) and 

p q 

(4.2) p d ( Z )= n (l+^na+M 2 " 

3=5+1 i=0 

for 9i = 27(2% + 1) + 2P+ 1 - 2?+ 1 ). Note that we take z = -KX, where K = 
&! + ••• + k m is the sum of all time steps in the sequence, rather than z — —KX, where 
K < K is the maximum time step in the sequence. To make the comparison fair, we 
take m = 5 for the Chebyshev polynomial and (p, q) = (3, 1) for the dyadic damping 
polynomial, which gives the same number of time steps (to = 2 q+1 — I + p — q = 5) 
for the two methods. 

From Figures [1] and [2J we see that Chebyshev damping can be slightly more 
efficient than dyadic damping; the cost is smaller for a given value of KXn (Figure [TJ 
and the stability interval on the negative real axis is larger (Figure [2]). However, 
with dyadic damping we don't need to assume that the eigenvalues of A in (13. 1| lie 
in a narrow strip along the negative real axis in the complex plane, as is needed 
with Chebyshev damping. (The stability region of the Chebyshev polynomial can be 
slightly modified to include a narrow strip along the negative real axis; see [11].) 
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Fig. 1. A comparison of the two polynomials for K\n = 64: we take m = 6 for the shifted 
Chebyshev polynomial P c (%) and (p,q) = (6,3) for the dyadic polynomial P^(x). With K = 1, the 
costs are 5.3 and 8, respectively. 
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Fig. 2. A comparison of stability regions for the two polynomials, with m = 5 for the shifted 
Chebyshev polynomial P c (z) (left) and (p,q) = (3,1) for the dyadic polynomial Pd{z) (right). 



5. The general nonlinear problem. We consider now the general nonlinear 
problem, 



(5.1) 



u(t) = f(u(t)) for t > 0, 
u(0) = u°, 
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where / : M. N x (0, T] — > M. N is a given bounded and differentiable function. The 
explicit Euler method for (|5.ip reads 

(5.2) U n = U n ~ l + k n f{U n ~ l ), 

where the function U = U(t) is piecewise constant and right-continuous with U n = 
U(t+). The exact solution u satisfies a similar relation, 

(5.3) u n = u n - x + f " /(«(*)) dt, 

with u 11 = u(t n ). Subtracting (|5.2[) and (|5.3[) , we find that the error e(t) — U(t) — u{t) 
satisfies 



e(tf)-e(Ci)=/ J(U(t)-u(t))dt= Jedt, 

where J is the Jacobian & of the right-hand side evaluated at a mean value of U and 
u. We conclude that the efficiency of the proposed method for the general nonlinear 
problem will be determined by the distribution of the eigenvalues of the Jacobian. 

6. Iterative methods. From another viewpoint, we may consider using an 
explicit-type iterative method for solving the discrete equations arising from an im- 
plicit method. The implicit cG(l) method with midpoint quadrature for the general 
nonlinear problem (|5.1|) reads 



3.1) U n = U n - X + k n f ' 



We can solve this system of nonlinear equations for U n using Newton's method, but 
the simplest and cheapest method is to apply fixed point iteration directly to (|6.ip : 
that is, 

'jjn-x , yni-x 

U nl = u n - L + k n f ' 



for I = 1, 2, . . . until convergence occurs with U n, ° = C/™ -1 . The fixed point iteration 
converges for a small enough time step k n , so the stability condition for a standard 
explicit method appears also in explicit-type iterative methods as a condition for 
convergence of the iterative solution of the implicit equations. 

To determine a stop criterion for the fixed point iteration, we measure the size of 
the discrete residual, 

r nl = ± (U nl -C/"- 1 )-/ / 

which should be zero for the true cG(l) approximation. We continue the iterations 
until the discrete residual is smaller than some tolerance tol > 0. Usually only a 
couple of iterations are needed. Estimating the error e(T) = U(T) — u(T) at final 
time T (see [8]), we have 

||e(T)|| < S(T) maxfc||i?|| +S°(T) max llrll, 

[0,T] [0,T] 
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where R(t) = U(t) — f(U(t)) is the continuous residual and S(T) and S°(T) are 
stability factors. For the test equation, we have S(T) < 1 and S°(T) < 1/A, which 
suggests that for a typical stiff problem we can take tol = TOL, where TOL is a 
tolerance for the error e(T) at final time. 
For the discrete residual, we have 



Thus, by measuring the size of the discrete residual, we obtain information about 
the stiff nature of the problem, in particular the eigenvalue of the current unstable 
eigenmode, which can be used in an adaptive algorithm targeted precisely at stabiliz- 
ing the current unstable eigenmode. We discuss this in more detail below in section [HJ 

7. Mult i- adaptive solvers. In a multi-adaptive solver, we use individual time 
steps for different components. An important part of the algorithm described in 
[8j |9] is the iterative fixed point solution of the discrete equations on time slabs. The 
simple strategy to take small damping steps to stabilize the system applies also in 
the multi-adaptive setting, where we may also target individual eigenmodes (if these 
are represented by different components) using individual damping steps. This will 
be explored further in another note. 

8. An adaptive algorithm. The question is now whether we can choose the 
time step sequence automatically in an adaptive algorithm. We approach this question 
in the setting of an implicit method combined with an explicit-type iterative solver as 
in section^ We give the algorithm in the case of the simple damping strategy outlined 
in section [2] with an extension to parabolic problems as described in section [3] 

A simple adaptive algorithm for the standard cG(l) method with iterative solution 
of the discrete equations reads as follows. 

1 . Determine a suitable initial time step ki , and solve the discrete equations for 
the solution U(t) on (£ , *i ) ■ 

2. Repeat on (t n _x,i n ) for n = 2, 3, . . . until t n > T: 

(a) Evaluate the continuous residual R n -i from the previous time interval. 

(b) Determine the new time step k n based on R n -\. 

(c) Solve the discrete equations on (i n _i,t n ) using fixed point iteration. 

In reality we want to control the global error, which means we also have to solve the 
dual problem, compute stability factors (or weights), evaluate an a posteriori error 
estimate, and possibly repeat the process until the error is below a given tolerance 
TOL > 0. The full algorithm is thus slightly more elaborate, but the basic algorithm 
presented here is the central part. See [1] for a discussion. 

We comment also on step 2(b): For the cG(l) method we would like to take 
k n = TOL/(5(T)||i? n _i||), but this introduces oscillations in the size of the time 




which gives 



(6.2) 




Jr 



.71,1-1 
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step. A small residual gives a large time step which results in a large residual, and so 
on. To avoid this, the simple step size selection has to be combined with a regulator 
of some kind; see [H [TU] or [5]. It turns out that a simple strategy that works well in 
many situations is to take k n as the geometric mean value 

(8.1) k n = - 2KK - 1 



where k n = TOL/(S(T)|| J R n _i||). 

Now, for a stiff problem, what may go wrong is step 2(c); if the time step k n is too 
large, the fixed point iterations will not converge. To be able to handle stiff problems 
using the technique discussed above, we propose the following simple modification of 
the adaptive algorithm. 

1 . Determine a suitable initial time step fei , and solve the discrete equations for 
the solution U(t) on (to,ii)- 

2. Repeat on (t n —i, t n ) for n = 2, 3, . . . until t n > T: 

(a) Evaluate the continuous residual Rn-i from the previous time interval. 

(b) Determine the new time step k n based on R n -\. 

(c) Solve the discrete equations on (t n -i,t n ) using fixed point iteration. 

(d) If 2(c) did not work, compute 

T _ 2 Hr'll 



k n Ik'- 1 !!' 

and take m = log(fc„L) explicit Euler steps with time step k — c/L and 
c close to 1. 

(e) Try again starting at 2(a) with n — > n + m. 

In the analysis of section[2j we had c = 1/2, but it is clear that the damping steps 
will be more efficient if we have c close to 1. An implementation of this algorithm in 
the form of a simple MATLAB code is available for inspection [7], including among 
others the test problems presented in the next section. 

We also note that, by (|8.1[) . we have k n < 2fc„_i, so, following the sequence 
of small stabilizing time steps, the size of the time step will be increased gradually, 
doubling the time step until we reach k n ~ K or the system becomes unstable again, 
whichever comes first. This automatically gives a sequence of time steps similar to 
that of dyadic damping described in section [3J with the difference being that most of 
the damping is made with the smallest time step. For a parabolic problem, we modify 
the algorithm to increase the time step more carefully, as determined by Table Q] with 
p given by p = log(fc„L)/ log(2). 

9. Examples. To illustrate the technique, we take a simple standard imple- 
mentation of the cG(l)-method (with an explicit fixed point solution of the discrete 
equations) and add a couple of lines to handle the stabilization of stiff problems. We 
try this code on a number of well-known stiff problems taken from the ODE literature 
and conclude that we are able to handle stiff problems with this explicit code. 

When referring to the cost a below, this also includes the number of fixed point 
iterations needed to compute the cG(l) solution on intervals where the iterations 
converge. This is compared to the cost ao for the standard cG(l) method in which we 
are forced to take small time steps all the time. (These small time steps are marked 
by dashed lines.) For all example problems below, we report both the cost a and the 
cost reduction factor aja^. 
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9.1. The test equation. The first problem we try is the test equation, 

u(t) + Xu(t) =0 for t > 0, 

(9.1) W n 
V ' u(0) = U°, 

on [0, 10], where we choose u = 1 and A = 1000. As is evident from Figure [3l the 
time step sequence is automatically chosen in agreement with the previous discussion. 
The cost is only a « 6 and the cost reduction factor is a/ao « 1/310. 

Note how the time steps are drastically decreased (step 2(d) in the adaptive 
algorithm) when the system needs to be stabilized and then gradually increased until 
stabilization again is needed. 

9.2. The test system. For the test system, 

ii(t) + Au(t) =0 for t > 0, 

(9.2) W W 

1 ' u(0) = u°, 

on [0, 10], we take A = diag(100, 1000) and u° — (1, 1). There are now two eigenmodes 
with large eigenvalues that need to be damped out. The dominant eigenvalue of A is 
A2 = 1000, and most of the damping steps are chosen to damp out this eigenmode, but 
some of the damping steps are chosen based on the second largest eigenvalue Ai = 100 
(see Figure |4]) . When to damp out which eigenmode is automatically decided by the 
adaptive algorithm; the bad eigenmode that needs to be damped out becomes visible 
in the iterative solution process. Since there is an additional eigenvalue, the cost is 
somewhat larger than for the scalar test problem, aw 18, which gives a cost reduction 
factor of a/ao ~ 1/104. 
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9.3. A linear nonnormal problem. The method behaves similarly even if we 
make the matrix A highly nonnormal. We now solve 



(9.3) 

on [0, 10], with 



u(t) + Au{t) 

u(0) = a 



for t > 0, 
o 



A = 



1000 




10000 
100 



and u = (1, 1). The cost is about the same as for the previous problem, a « 17, but 
the cost reduction factor is better: a/ao ~ 1/180. The solution and the time step 
sequence are given in Figure [5j 

9.4. The HIRES problem. The so-called HIRES problem (High Irradiance 
RESponse) originates from plant physiology and is taken from the test set of ODE 
problems compiled by Lioen and de Swart [5J. The problem consists of the following 
eight equations: 

111 = -1.71ui + 0.43u 2 + 8. 32u 3 + 0.0007, 

112 = 1.71«i — 8.75i*2, 
u 3 = -10.03u 3 + 0.43u 4 + 0.035m 5 , 
in = 8.32u 2 + 1.71u 3 - 1.12u 4 , 
it 5 = -1.745u5 + 0.43u6 + 0.43u 7 , 
u 6 = -280.0u 6 -it 8 + 0.69u 4 + 1.71u5-0.43u 6 + 0.69w7, 

1*7 = 280.01*6^8 ~ I.8IM7, 

iis = -280.0u 6 u 8 + 1.81u 7 , 



(9.4) 
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Fig. 5. Solution and time step sequence for l|9.3ll . ct/cto £S 1/180. 

together with the initial condition u° = (1.0, 0, 0, 0, 0, 0, 0, 0.0057). We integrate over 
[0, 321.8122] (as specified in [6]) and present the solution and the time step sequence 
in Figure [6l The cost is now a ~ 8, and the cost reduction factor is a/ao ~ 1/33. 

9.5. The Akzo Nobel problem. The next problem is a version of the "Chem- 
ical Akzo-Nobel" problem taken from the ODE test set [5] , consisting of the following 
six equations: 

id = -2rr + r 2 - r 3 - r 4 , 
ii2 — — 0.5ri — r4 — 0.5r5 + F, 
(9.5) " 3 = ^-^+^3, 

ii5 = r 2 -r 3 +r 5l 
, ue = -rs, 

where F — 3.3- (0. 9/737— u^) and the reaction rates are given by n = 18.7-uf y 7 ?^, ?*2 = 
0.58-M3W4, f3 = 0.58/34.4 -ui«5, r 4 = 0.09-mim|, and r§ = 0.42 • Mgy 7 ?!^. We integrate 
over the interval [0,180] with initial condition u° = (0.437,0.00123,0,0,0,0.367). 
Allowing a maximum time step of fc ma x = 1 (chosen arbitrarily), the cost is only 
a ~ 2, and the cost reduction factor is a/ao ~ 1/9- The actual gain in a specific 
situation is determined by the quotient of the large time steps and the small damping 
time steps, as well as the number of small damping steps that are needed. In this 
case, the number of small damping steps is small, but the large time steps are not 
very large compared to the small damping steps. The gain is thus determined both by 
the stiff nature of the problem and the tolerance (or the size of the maximum allowed 
time step). 
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9.6. Van der Pol's equation. A stiff problem discussed in the book by Hairer 
and Wanner is Van der Pol's equation, 



u + ii{u 2 — l)u + u = 0, 
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which we write as 



(9.6) 



We take /i = 1000 and compute the solution on the interval [0, 10] with initial condition 
u a = (2,0). The time step sequence behaves as desired with only a small portion of 
the time interval spent on taking small damping steps (see Figure [8]). The cost is now 
a « 140, and the cost reduction factor is a/ao ~ 1/75. 

9.7. The heat equation. A special stiff problem is the one-dimensional heat 
equation, 

ii(x,t) - u"(x,t) = f(x,t), 16(0,1), t>0, 
u(0,t) = u(l,t) = 0, t>0, 

u(x,Q) = 0, a; 6 [0,1], 

where we choose f(x,t) — f(x) as an approximation of the Dirac delta function at 
x = 0.5. Discretizing in space, we obtain the ODE 

[ ' ' «(o) = o, 

where A is the stiffness matrix. With a spatial resolution of h = 0.01, the eigenvalues 
of A are distributed in the interval [0,4 • 10 4 ] (see Figure [9]). 

Using the dyadic damping technique described in section [3] for this parabolic 
problem, the cost is a sa 2000, with a cost reduction factor of a/ao w 1/31. 
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10 20 30 40 50 60 70 80 



Fig. 9. Solution and time step sequence for l|9.7jl , a/ao 1/31. Note that the time step 
sequence in this example is chosen slightly differently than in the previous examples, using the 
technique of dyadic damping discussed in section [3] This is evident upon close inspection of the 
time step sequence. 
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20 40 60 80 100 




95 96 97 98 99 100 




Fig. 10. Solution and time step sequence for l|9.8j l. a/ao = 1. 



9.8. A nonstiff problem. To show that the method works equally well for 
nonstiff problems, we finally consider the following simple system: 



(9.8) 



ill = 5u 2 , 
u 2 = 



With the initial condition u° = (0, 1), the solution is u (t) = (v / 5sin( v / 50:Cos(v / 5)). 
Since this problem is nonstiff (for reasonable tolerances), no stabilization is needed, 
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and so the solver works as a standard nonstiff cG(l) solver with no overhead. This 
is also evident from the time step sequence (see Figure 1 10[) which is chosen only to 
match the size of the (continuous) residual. 
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